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Abstract 

o 

The hydrodynamic state of an impurity immersed in a low density granular gas is analyzed. 
Explicit expressions for the temperature and density fields of the impurity in terms of the hydro- 
^ , dynamic fields of the gas are derived. It is shown that the ratio between the temperatures of the two 

components, measuring the departure from energy equipartition, only depends on the mechanical 

properties of the particles, being therefore constant in the bulk of the system. This ratio plays an 

- 

important role in determining the density profile of the intruder and its position with respect to the 

o. 

gas, since it determines the sign of the pressure diffusion coefficient. The theoretical predictions are 

compared with molecular dynamics simulation results for the particular case of the steady state of 

. an open vibrated granular system in absence of macroscopic fluxes, and a satisfactory agreement 

O . 

£NJ , is found. 
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I. INTRODUCTION 



As a consequence of the inelasticity of collisions, granular fluids exhibit a series of behav- 
iors that are in sharp contrast with those of molecular fluids. One of them is the absence of 
energy equipartition in a granular mixture of mechanically different species. The granular 
temperatures of the components of the mixture, defined from the average kinetic energies, 
are different. Although this possibility was already pointed out many years ago it has 
not been until recently that a systematic study of the effect has been started. For the homo- 
geneous state of a freely cooling binary mixture of inelastic spheres, an explicit expression 
for the ratio of temperatures of the two components has been obtained from an approxi- 
mate solution of the kinetic Enskog equations 2]. The accuracy of this prediction for weak 
dissipation and low density has been confirmed by Molecular Dynamics (MD) simulations 

m 

The above homogeneous cooling state (HCS) is not accesible experimentally. In order 
to maintain a granular system fluidized, an external energy supply is required. This is 
often carried out by means of vibrating walls or external fields, that generate macroscopic 
gradients in the system because of inelasticity. At a theoretical level, homogeneously driven 
granular fluids have also been considered. In these models, stochastic forces injecting energy 
are added, then allowing a system otherwise isolated to reach a steady state. The lack of 
energy equipartition in the steady state of homogeneously driven granular mixtures has been 
also analyzed J], by extending the methods of Q]. Nevertheless, the relationship between 
this kind of driving and actual experiments is uncertain. The non-equipartition has also 
been confirmed by MD simulations of simple shear flows 5] and of vibrated granular gases 
(in absence of gravity) 

Experimental evidence of the coexistence of different temperatures in strongly vibrated 

n n 

granular mixtures has been reported in both two-dimensional p]] and three-dimensional [8( 
systems. Moreover, the dependence of the temperature ratio on the several parameters 
characterizing the system was investigated. Quite surprisingly, it has been observed that 
the above ratio remains practically constat in the bulk of the system, in spite of the rather 
large gradients exhibited by each of the partial temperatures. 

The lack of energy equipartition has consequences on the expressions of the transport 
coefficients of a granular mixture j^]. Therefore, it must also affect the shape of the hydro- 
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dynamic profiles and, in particular, the density distribution of each of the components in 
inhomogeneous mixtures. Consequently, it seems clear that the temperature difference can 
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play a role in the segregation phenomena occurring in vibrated granular mixtures 

Here, the particular case of a low density binary mixture in the tracer limit, i.e. when 
the mole fraction of one of the components is very low, will be considered. The simplicity of 
the system allows a detailed and controlled discussion of the hydrodynamic profiles of the 
tracer component, that can be expressed analytically in terms of those of the other (excess) 
component. The starting point are the kinetic Boltzmann equations for the mixture, and 
the analysis is based on the Chapman- Enskog method. The theory is first formulated for a 
general kind of states and later on particularized for the steady state of an open vibrated 
system. Then, the theoretical predictions are compared with MD simulation results, and a 
quite satisfactory agreement is found over a wide range of values of the parameters defining 
the mechanical properties of the system. A short summary of some of the results presented 
here was given in ref. 1^ . 

The plan of the paper is as follows. In Sec. m some on the fundamental ideas on which 
the Chapmann- Enskog method, as applied to a one-component system, is based are shortly 
reviewed. They are used in Sec. IHIl where an explicit expression for the ratio between the 
local temperatures of the two components is derived. This ratio only depends on the me- 
chanical properties of the particles, being independent of the gradients of the hydrodynamic 
fields. An extensive comparison with MD results for a vibrated granular system in presence 
of gravity is carried out. The shape of the density profile for the tracer component is inves- 
tigated in Sec. IIV| where an expression for the mass flux relative to the local flow is derived 
to first order in the gradients of the mixture fields (Navier-Stokes order). For the particular 
case of the steady state of an open vibrated system, an expression for the density profile 
of the tracer component in terms of the fields of the excess component is obtained. Since 
analytical forms for these fields are known, the expression can be integrated numerically. 
Some details of this are given in the Appendix. Again, the theoretical predictions compare 
well with MD simulation results, although some significant discrepancies show up as the 
difference in mass of the particles is increased. Special emphasis is put on the relationship 
between the breakdown of energy equipartition and the relative position of the centers of 
mass of the two species. Finally, Sec. |V| contains some further discussion of the results, as 
well as some comments on their applicability to arbitrary states. 
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II. CHAPMAN-ENSKOG DESCRIPTION OF A GRANULAR GAS 



As already said, the aim of this work is to analyze some properties of the hydrodynamic 
profiles of a dilute granular binary mixture in the tracer limit, i.e. when the mole fraction 
of one of the species is very small. More precisely, the emphasis is put on the relationship 
between the profiles of the tracer component and those of the mixture in an arbitrary 
state. It is assumed that the concentration of the tracer component is so small that its 
presence does not affect the state of the granular mixture that is, therefore, determined by 
the state of the other (excess) component. At a kinetic theory level, this implies that the 
evolution of the one-particle distribution function of the excess component obeys a closed 
nonlinear Boltzmann equation. Moreover, in the evolution of the distribution function of 
the tracer component, the mutual interactions between their particles can be neglected, as 
compared with the interactions with the particles of the excess component. Consequently, 
its distribution obeys a linear Boltzmann-Lorentz equation. This is formally equivalent to 
consider an impurity or intruder immersed in a dilute granular gas, and this will be the 
terminology used in the following. Let us start by shortly reviewing some basic aspects of 
the Chapman-Enskog method applied to a one-component granular gas that will be needed 
in the following. 

Consider a low density gas composed by smooth inelastic hard spheres (d = 3) or disks 
(d = 2) of mass m and diameter a. The inelasticity of collisions is modelled by a velocity- 
independent coefficient of normal restitution a, defined in the interval < a < 1. The 
particles are submitted to an external force of the gravitational type, so the force acting on 
each particle has the form F = —mg e z , where g is a constant and e z the unit vector in the 
positive direction of the z axis. It is assumed that in the low density limit, the time evolution 
of the one-particle distribution function of the gas, f(r,v,t), is accurately described by the 
nonlinear inelastic Boltzmann equation 



where J is the Boltzmann collision operator describing the scattering of pairs of particles, 




(1) 




xg-a[a 2 f(r,v' ,t)f(r,v' 1 ,t) 
-f(r,v,t)f(r,v 1 ,t)\. 



(2) 
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Here is the Heaviside step function, <r is the unit vector pointing from the center of 
particle 1 to the center of the other colliding particle at contact, and g = v — v\ the relative 
velocity. Moreover, v' and v[ are the precollisional or restituting velocities, i.e. the initial 
values of the velocities leading to v and V\ following the binary collision defined by <r, 

v = v--(g-v)*, 

v[ = Vl + ^-^{g-a)a. (3) 

Macroscopic fields number density, n(r,t), flow velocity u(r,t), and granular temperature, 
T(r,t), are defined in the usual way as velocity moments of the distribution function, 



n{r,t) = J dvf(r,v,t), 
n(r,t)u(r,t) = J dvvf(r,v,t), 



/TflV 2 
dv—f(r,v,t), (4) 

where V = v — u is the peculiar velocity. Balance equations for the above fields are directly 
obtained from the Boltzmann equation, 

d t n + V • (nu) = 0, (5) 

d t u + u ■ Vn+ (mn) _1 V • P + g e z = 0, (6) 
d t T + u-VT + 2{nd)-\W ■ q + P : Vu) + T( = 0. (7) 
The pressure tensor P and the heat flux q are given by 

P = J dvmVVf(r,v,t), (8) 

q = Jdv^V 2 Vf(r,v,t), (9) 

respectively. The cooling rate, (, in Eq. (jZJ is due to the energy dissipation in collisions and 
it is also a functional of the distribution function, 

C(r,t) = -(nTd)- 1 J dvmv 2 J[v\f,f}. (10) 

The balance equations ©-(0) only become a set closed hydrodynamic equations for the 
fields once P, q and ( are expressed as functionals of them. This requires to find a solution 
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to the Boltzmann equation such that all the space and time dependence occurs through n, u, 
and T. In the Chapman-Enskog procedure, this normal solution is generated by expressing 
it as a series expansion in a formal non- uniformity parameter e, 

/ = /(°) + e/« + 6 2 / (2) + ..., (11) 

where each factor of e means an implicit gradient of a macroscopic field. Use of the above 
expansion in the definitions of the fluxes and the cooling rate gives a corresponding expansion 
for them. Finally, the time derivatives of the fields are also expanded in the gradients, 
d t = + edf + . . ., by means of the balance equations. Then, it is 

$°>n = 0, di 0) u = 0, di 0) T = -TC {0 \ (12) 

where the zeroth order cooling rate is given by 

C (0) (r,t) = -{nTd)- 1 J dvmV 2 J[v\f^\f Q) }. (13) 

In this way, to zeroth order in e it is obtained: 

-C (0) r^ = J[H/ (0) ,/ (0) ]. (M) 

The solution to this equation is chosen such that it verifies the conditions 

J dvvf(°\r,v,t) = n(r,t)u(r,t), 

J dv^f(°\r, v, t) = n(r, t)T(r, t), (15) 

i.e., it leads to the same hydrodynamic fields as the complete distribution function. This 
choice is consistent with the solubility conditions for the set of equations generated by the 
Chapman-Enskog method Q. Thn s /<»> is easily generated from the distribution function 
for the homogeneous cooling state (HCS) of a dilute gas 14| by substituting the homogeneous 
density and temperature fields by the local values n{r,t) and T(r,t), and replacing the 



velocity v by the peculiar one, V(r,t). Therefore, it has the form 14\ 

f(°\r,v,t) = nvT d (t) X (V/v T ), (16) 
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where v%(r,t) = 2T(r,t)/ m is a local thermal velocity. As a consequence, Eq. ()14|) can be 
rewritten as 

^^•(^ (0) ) = J N/ (0) J (0) ]- (17) 

Using Eq. (|TH|l it is easily obtained: 

d^n = -V • (nu), (18) 

9j u = — u ■ Vu — (mn) _1 Vp — go^, (19) 

d[ 1] T = -u-VT- • u. (20) 

In the last expression, use has been made of the fact that = because of symmetry 
— o ns Q. Moreover,^™ pressure of the gas. The procedure can be 
continued and Navier-Stokes hydrodynamic equations for a dilute granular gas, with explicit 
expressions for the transport coefficients in the first Sonine approximation, have been derived 



expressi 

HQ. 



III. TEMPERATURE OF THE IMPURITY 

Let us suppose now that an impurity or intruder of mass mo and diameter <To is added 
to the system. It is assumed that the presence of the intruder does not affect the state of 
the gas, so that the one-particle distribution function of the gas particles is still determined 
by the nonlinear Boltzmann equation Moreover, the macroscopic flow velocity and 

temperature for the mixture formed by the gas plus the intruder are the same as those for 
the gas alone, i.e. they are given by Eqs. (jlj). 

The distribution function of the intruder, fo(r, v, t), obeys the linear Boltzmann-Lorentz 
equation 

(dt + v-V- g -^pj f (r, v, t) = J [v\f 0l /], (21) 
where the collision operator now is 



a^f (r,v",t)f(r,v>;,t) 



xg ■ a 

-f (r,v,t)f(r, Vl ,t)}. (22) 
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Here, a = (a + <Jo)/2 and «o is the coefficient of normal restitution for impurity-gas collisions. 
The precollisional velocities in this case are given by 

„ (l + «o)A, 
v = v {go-jo-, 

(l + a )(l- A) . 
v 'l = Vl + S °A L( g . (23) 

a 

with A = mj (m + m ). The number density for the intruder is 

n (r,t) = J dvf (r,v,t). (24) 

The Chapman- Enskog procedure when applied to a mixture, assumes the existence of a 
normal solution of the Boltzmann equations for the mixture in which all the space and time 
dependence of the one particle distribution function of each of the s pec ies occurs through 
a functional dependence on the hydrodynamic fields of the mixture [13|. These fields can 
be chosen in different, equivalent ways. Here they are taken to be the concentration of 
the impurity xq = n Q /(n + n) — > n /n, the pressure p, the local flow velocity u, and the 
temperature T, i.e. 

f (r,v,t) = f [v\x (r,t),p(r,t),u(r,t),T(r,t)]. (25) 

Then, an expansion similar to the one given in Eq. (fTTj) is considered for the distribution 

function of the impurity, f = + e/o + t 2 f^ + To lowest order in e, Eq. (|2Tj) 

becomes 



-C ( °'(p| + rA)# = JoN/„ <0) , /<»>]. 

Taking into account that dimensional analysis requires that /q is of the form 



(26) 



/o (0) =*of vj d Xo (V/vr), (27) 



Eq. (j26|) is seen to be equivalent to 

2 dV 

The solution of this equation is chosen such that 



(v/ (0) )=Jo[^|/o (0) ,/ (0) ]. (28) 



J dvf(°\r,v,t) = J dvf (r,v,t)=no(r,t), (29) 
J dv vti°\r,v,t) = n (r,t)u(r,t), (30) 



J dvmV 2 fo°\r,v,t) = J dvmV 2 f (r , v , t) 

= dn (r,t)T (r,t). (31) 

The last equality in Eq. (|3*TJl defines the local temperature, T , of the intruder. The above 
requirements are consistent with the solubility conditions of the equations generated by the 
Chapman-Enskog method. It is worth to stress that Eq. (|2*Hj) holds with independence of 
the specific form of the hydrodynamic fields of the mixture. Thus multiplication of that 
equation by moV 2 and integration over V yields 

C (0) (r,t) = Ci 0) M), (32) 

where 

ti°\r,t) = -KTodr 1 / dVm V 2 J [v\f!?\fW} (33) 

is the lowest order in the gradients of the cooling rate for the temperature of the intruder. 

The evaluation of the cooling rates and Q requires to solve Eqs. (j!4j) and (|28j) . 
This can be done in a systematic way by expanding /(°) and /q m terms of an ensemble 
of orthogonal polynomials 2!]. Here we will consider a leading order approximation that is 
expected to give quite accurate results, at least for not very strong inelasticity. The zeroth 
order distributions are approximated by Gaussians, 



/2vrT\ 



-d/2 



f^(r,v,t)=n(— e"^, (34) 
V m / 

ti°\r,v,t)=nJ 2 -^Y d/2 e--^. (35) 



V m 

Note that these expressions are consistent with the conditions (fPfj) and ([29 )1 - (|31 |) . Employing 
them, it is straightforward to calculate the cooling rates. The technical details needed to 



evaluate the integrals have already been discussed several times 
be reproduced here. The result is 



7j , and they will not 



C (0) M) V2^l 2 



n°>* = 2 = U^l (i _ a 2 ) (36) 



^ - 7 C# m = W + ( 1 - h l -±±) , (37) 



where h = m(l + ao)/2(m + m Q ), Vq is a dimensionless collision rate, 

8h7r (d-l)/2 /grxd-l 



T(d/2)d \a 



(38) 



and is the ratio of the mean square velocities for the intruder and fluid particles, 

= mT (z) 
m Q T(z)' 

Substitution of Eqs. (J36|) and pTjl into Eq. (J32|) gives 



(39) 



where the parameter /3 is given by 

" s vr(?) ■ (41) 

The solution of this equation provides the expression for T Q (r,t)/T(r,t). There is a unique 
real solution for all allowed values of h and f3. For elastic collisions (a = = 1), the 
solution is <j) = m/rriQ and energy equipartition follows. If only the intruder-gas collisions 
are inelastic (a = 1), the result = h/(l — h) derived by Martin and Piasecki for an inelastic 
impurity in an equilibrium elastic gas is recovered 

The simplicity of Eq. (|40|) may appear as a surprise, specially taking into account its 
generality. No particular state of either the gas or the intruder has been considered and, 
nevertheless, the temperature ratio is given as a function of only the mechanical properties 
of the particles. The existence of an expression for the temperature ratio is a consequence of 
the assumption that there is a hydrodynamic description of the mixture in terms of only xq, 
p, u, and T. Therefore, although the partial temperatures of the intruder and the gas are 
different, they are not needed to specify the macroscopic state of the mixture. Nevertheless, 
the theory provides a relationship to determine them. The fact that Eq. ()4U|) does not 
involve any of the other hydrodynamic fields follows from the tracer and dilute limits we are 
considering. 

The above relation between the partial temperatures is not modified when applied to 
the HCS, consistently with previous results This does not imply by itself any kind 

of local HCS approximation, in the same way as energy equipartition in molecular systems 
does not imply local equilibrium, and it is valid independently of the gradients present in 
the system. On the other hand, the above results do not hold if the characterization of the 
macroscopic state of the mixture also requires to specify the partial temperatures of the gas 
and the intruder and, therefore, the set of independent hydrodynamic variables (and also 
the number of equations) must be expanded. This seems to be the kind of approximation 
followed in 
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In order to check the accuracy of the theoretical prediction given by Eq. (|40|). we have 
performed MD simulations with an event driven algorithm. The particular situation we have 
considered is a two-dimensional system confined in a rectangular box of width W, open on 
the top, and submitted to an external field of the gravitational type, as defined above Eq. 

Periodic boundary conditions are enforced in the direction perpendicular to the field. 
To maintain the system fluidized, energy is continuously supplied through the wall located 
at z — 0, which is vibrating. For simplicity, and given that we are interested in the bulk 
properties of the system, the wall moves in a sawtooth way, with very small amplitude and 
high frequency, so that it can be considered that all the particles colliding with the wall find 
it at z = and with the same velocity v w = vy^e z 2(3, 21]. The particle- wall collisions are 
elastic. 

Under the above conditions, it has been proven that a one-component granular gas reaches 
a stationary state with only gradients in the ^-direction j^]. Nevertheless, this state becomes 
unstable and develops transversal inhomogeneities when the size of the vibrating wall is 



larger than a critical value 



22|. In fact, the same happens for a closed system in absence of 



gravity 23, 24]. In all the simulations to be reported in the following, the value of W has 
been taken small enough so the system stays homogeneous in the transversal direction. 

The above steady state is highly inhomogeneous. An example is given in Fig^ where 
the density and temperature profiles for a monodisperse system of N = 359 particles with 
a = 0.95 are shown. In order to plot the temperature profile on the same scale as the density 
one, the former has been scaled with an arbitrary value Ta- Here and in the following, the 
units defined by m — 1, a — 1, and go = 1 will be used. In these units, the data in the 
figure correspond to W = 50 and vw = 5. These values lead to a fluidized state without 
transversal inhomogeneities, and with a density low enough as to expect the Boltzmann 
equation to give an accurate description of the system. It is seen that the density profile 
exhibits a maximum at which n ~ 0.075. The existence of a density maximum is a general 
: eature of open vibrated systems as long as the number of monolayers at rest is large enough 
On the other hand, the temperature profile presents a minimum which is related with 
the existence of a term proportional to the density gradient in the expression of the heat 



flux 



21 



3. 

Figures |2H show the profile of the temperature ratio T (z)/T(z) for different values of the 
mechanical properties of the particles. It is observed that, aside statistical fluctuations, the 
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o n 

aT/T a 




FIG. 1: Density and temperature vertical profiles for a vibrated granular gas in the steady state 
considered in the text. The number of particles is N = 359, the restitution coefficient a = 0.95, 
and the size of the vibrating wall W = 50a. The symbols are from MD simulations and the solid 
lines the fit to the hydrodynamic prediction discussed in the Appendix. 
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FIG. 2: Temperature ratio profiles Tq(z)/T{z) for several values of the restitution coefficient ao, as 
indicated in the figure. The values of the other parameters are a = 0.95, ctq = a, and mo = m/2. 
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FIG. 3: Temperature ratio profiles Tq(z)/T(z) for several values of the mass ratio mo/m, as 
indicated in the figure. The values of the other parameters are a = 0.95, ao = 0.7, and o"o = o. 



ratio remains practically constant in a wide region of the system, even if T (and T ) varies 
significantly. In Fig. El the influence of the inelasticity of the gas-intruder collisions on the 
temperature ratio is analyzed. As expected, the deviation from equipartition increases as ao 
decreases, becoming more different from a. The influence of the mass ratio m Q /m is studied 
in Fig. El where it is seen that the temperature ratio significantly deviates from unity when 
itlq = m but ao differs significantly from a. This is in contrast with previous findings in 
both experiments tJ and MD simulations of vibrated granular mixtures (not in the tracer 
limit), probably because the several coefficients of normal restitution in those works were 
too close one another as to observe this deviation. The figure also shows that the boundary 
layer near the vibrating wall, defined as the region in which the hydrodynamic prediction is 
not verified, increases as the ratio mo/m increases. This can be explained as a consequence 
of the larger increase of the kinetic energy of the impurity than of the gas particles when 
colliding with the vibrating wall. In addition, the collisions of the intruder with the gas 
become less efficient in "thermalizing" the former. The same explanation applies to the 
strong fluctuations at large heights. Finally, the sensitivity of Tq/T on the diameter ratio is 
illustrated in Fig. HI from where it can be concluded that it is quite weak. 

A quantitative comparison between the theoretical prediction given by Eq. ()4())1 and the 
MD results is presented in Figs. ElandlHl According with Eq. pUj) . for a given value of /?, 



13 



V 

1.5 



0.5 



o o a=2a 

H~~Q O =0.1O 



BB /■ 
B / 'rB 




150 



FIG. 4: Temperature ratio profiles Tq(z)/T(z) for several values of the diameter ratio <7o/cr, as 
indicated in the figure. The values of the other parameters are a = 0.95, «o = 0.6, and mo = m. 

the value of is a function of only the parameter h. Then, the comparison will be presented 
by considering series of simulations in which the value of the coefficient of restitution of the 
gas, a, and the diameter ratio, <j/<Jo, and therefore ft are kept constant, while the mass ratio 
m /m and the restitution coefficient and, therefore h are changed. In Fig. two series of 
data, one corresponding to cr = a, a = 0.8 (ft = 0.0636) and <j = a, a = 0.95 (ft = 0.0172) 
are presented. In each case, several values of the mass ratio have been considered, as 
indicated in the figure. In addition, the value of the coefficient of restitution for gas-intruder 
collisions has also been varied in the range 0.2 < cto < 0.99. To keep the figure readable, 
we have used the same symbols for the data corresponding to the same mass ratio although 
different values of q;q. A quite good agreement between theory and simulation is observed, 
specially for the largest value of a. The small but systematic discrepancies for a = 0.8 and 
mo/m = 3 are probably due to the failure of the hydro dyna mic Navier-Stokes equations to 



2l|. 



accurately describe the state of the g clS, clS discussed in 

The validity of the theory to describe the dependence on the diameter ratio is checked 
in Fig. |BJ Again, the same symbol is used for the MD results corresponding to situations 
differing only in the value of the coefficient of restitution a . As before, a weak dependence 
of the temperature ratio on the diameter of the particles is observed. It can be concluded 
that Eq. (|4"U|) provides an accurate description of the temperature ratio over a wide range 
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FIG. 5: Ratio of the mean square velocities eft as a function of the dimensionless parameter h defined 
in the main text. In all cases a = <jq. The solid and dashed lines are the theoretical prediction for 
a = 0.95 (/? ~ 0.0172) and a = 0.8 {(3 ~ 0.0636), respectively, while the empty and filled symbols 
are from the MD simulations for the same two values of a. The same symbol is used for simulation 
results differing only in the value of the restitution coefficient ceo- 

of values of the mechanical parameters characterizing both the gas and the intruder. 




IV. DENSITY PROFILE FOR THE IMPURITY 

An equation for the number density of the impurity no is directly obtained from Eq. 1)21 J). 

d t n + V ■ (no«) + rriQ 1 V • j = 0, (42) 
where jo{r,t) is the mass flux for the impurity relative to the local flow u(r,t), 

io(r, t)=m J dv Vf Q (r, v, t). (43) 



The isotropy of the function xo implies that the zeroth order of the flux mass for the 
aisf 

determine /q . To first order in e, Eq. (12 lj) reads 



intruder vanishes, i.e. = 0. Then, we proceed to study the next order, that requires to 



cfV (1) - ^o[/o (0) ,/ (1) ]-Jo[/ (1) ,/ {0) ] 



-BPfF-v-VfP+ga-g-fS*. (44) 
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FIG. 6: Ratio of the mean square velocities <p as a function of the dimensionless parameter h 
defined in the main text. In all cases a = 0.95. The solid and dashed lines are the theoretical 
prediction for ao/a = 2 {(3 ~ 0.0115) and &o/a = 0.1 (fi ~ 0.0313), respectively, while the empty 
and filled symbols are from the MD simulations for those two values of the diameter ratio. The 
same symbol is used for simulation results differing only in the value of the restitution coefficient 

The action of the time derivatives and on the hydrodynamic fields is given by Eqs. 
(H2J) and ([T8 |) -(j20 j) . respectively. Then, evaluating the right hand side of Eq. gives 



where the coefficients of the gradients of the fields on the right hand side are given by 

AO) 



^ (0) /o (1) -Jo[/o (0) ,/ (1) ]-^o[/o (1) ,/ (0) ] 




(45) 




(46) 




(47) 




(48) 




(49) 



16 



Here, I is the unit tensor in d dimensions and xq was defined above Eq. (|25j). Note that the 
external field does not occur in Eq. (|44|). This is because of the particular form for it we 
are considering (gravitational type). To proceed, we need the expression of It has been 
derived in ref. |l5| and reads 



/« = B ■ Vp - C ■ VT + V : Vu, 



(50) 



where the coefficients B, C, and D are functions of the peculiar velocity V and the hydro- 
dynamic fields. Then, the solution of Eq. (J45j) is of the form 



/ (1) = Aq ■ Vx + B ■ Vp + C ■ VT + V : Vu. 



(51) 



Substitution of Eqs. (joUj) and (joTj) into Eq. (|4*5|l and identifying coefficients of independent 
gradients yields 



c 



(0) 



9 MT d \ A 

y OP OP I 



- MA , /<°>] = A, 



(52) 



-Co 



' d m d \ 

"dp- + T of) +2 



Bo-Mf^^j-MBoJ 



(0)1 



B + -^C , 
P 



(53) 



-Co 



' d 3 \ 1 



Co- J [/o (1) ,C]- J [Co,/ (0) ] 



C - -^B , 



- Jo[/«S 0) ,^]- Jo[^o,/ (0) ] 
D(V). 



(54) 



(55) 



Use of Eq. (joT|) into Eq. (^H)) taking into account symmetry considerations, gives the ex- 
pression for the mass flux of the intruder to first order in the gradients, 



Tfl TflTl 

j = -m DVx - -D p Vp - —DVT, 



where 



is the diffusion coefficient, 



D 



D = -- J dvV ■ Aq 

Tm 



md 



dvV ■ Bo 



(56) 



(57) 



(58) 
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is the pressure diffusion coefficient, and 



D' 



Tm 
mnd 



dvV-Cn 



(59) 



is the thermal diffusion coefficient. 

In the following, the quantities Aq, Bo, and Co will be evaluated in the first Sonine 
approximation, i.e. we consider 



U \ 

Bo 
\C J 



/o ( V 



fa \ 

bo 
\c J 



(60) 



/ 

V 



f(°)V 



b 



(61) 



with /q°' ) being the Maxwellian given by Eq. (|35|). Consistently, B and C are approximated 
by 

' B 
C 

and it happens that both b and c are zero (see Appendix C in ref. |l5j]). Substitution of 
the above expressions into Eqs. (p^ - lpy^) gives a set of closed equations for a , bo, and Co- 
Multiplication of these equations by m V and integration over V yields 

n T 



-C (0) P~ + T— 



+ v 



c (0) [P^ + T 

\ op 







dT 



a n To 



2C (0) + v 



b n T 



(62) 



n o T o (0-l) T(Wcon T 

I 



V 



(63) 











dp 



K (0) MqT 
2T 



v 



c n T 



where 



v 



mo 



j >iw .. ioV.it -r" 



(64) 



(65) 



n T d 

is a collision frequency. It can be easily evaluated by using standard integration techniques 
with the result 

z/= ^ e(1 + ao) (l-A) 1/2 (l + 0) 1/2 (66) 
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where 



-<7 71 



m 



A V2 



(67) 



r(d/2)d 

is the elastic limit. An expression for v for the case d — 3, valid for arbitrary concentrations 
of the two components of the mixture and incorporating some non-gaussian contributions, 
has been obtained in p. Eq. (|6T)|) is consistent with that result. 

^From dimensional analysis, it follows that a n T ^ T 1 / 2 , b n T oc T 1 / 2 /p, and c n T cx 
ji-i/2 Then, the temperature and pressure derivatives in Eqs. (jH^-tEHI) can be evaluated 
and the equations become 



v 



3 £(0) \ 



feo 



a 



o > 



+ 



(o) 



V 



Co, 



c 



K (o) 



2Tv 



fen. 



Therefore, 



a = - 



fen = ~ 



SC V - 



' 3W0) a(o) 

V ; h 



Co 



2 

K (0) 

2TV 



2// 



feo- 



(68) 
(69) 
(70) 

(71) 
(72) 
(73) 



The expressions of the transport coefficients follow by using the above results into Eqs. 



D r , 



D 

n T ( 



nT 
-1) 



v 



C (oy 



-i 



3£(0) £(0)2 



mncj) 
D' 



2v 



2v p - 



for d 



(74) 
(75) 
(76) 

3 and without 



These expressions are consistent with those reported in ref. 
external field. 

Let us now particularize the above results for the steady state considered in the previous 
Section, i. e. an open vibrated system with no macroscopic flow and gradients only in the 
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z direction. In this case, Eq. ()42|) implies that j z = constant = 0, since the intruder mass 
flux must vanish for z — > oo. Thus Eqs. (j^HJ) and (f71 )l -f7S jl yield 

^ 2z/ J dz V 2i/ 2^ 2 j 

(dlnp C (0) 51nT\ 

Using Eqs. (jHSj) and (ffifij) . and also the equation determining 0, Eq. (flUj) . it is found 



( (0) 



1 



(78) 



This shows that C / v is a function of only /3 and /i, i.e. of the mechanical properties 
of the particles (0 is given in terms of them by Eq. f40|) ). but it does not depend on the 
hydrodynamic fields. Therefore, Eq. (fTTj) can be expressed in compact form as 

<9mxo 



dz 



-X(h : (3)g(z), (79) 



where 



and 



<91np C (0) <91nT 

^ ) = ^T-^7^- (81) 

Integration of Eq. (|T9"j) directly provides an expression for the probability density Pq{z) of 
finding the impurity at height z, 

P Q {z) = Cn(z) exp (-A J* dz'g(z')j , (82) 

with C being the normalization constant. 

The hydrodynamic profiles for the steady state of an open vibrated dilute gas (without 



impurity) where studied in [2jj and the results are shorty summarized in the Appendix. 
There, it is also indicated how the numerical evaluation of the integral appearing on the 
right hand side of Eq. (|82|) is implemented. 

In Figs. EKlthe theoretical and simulation results for Nxq and Pq{z) are compared for 
different values of the parameters. In all cases, those of the gas are the same as in Fig. 
i.e. a = 0.95, iV = 359, and W = 50, and the diameter ratio is (To/a = 1. On the other 
hand, in Fig. [7|it is m /m = 1/2, in Fig. |H1 m /m = 1, and in Fig. El m /m = 3. Finally, 
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100 200 300 

z/o 



FIG. 7: Normalized molar fraction Nxq (upper plot) and density profile Po(z) of the impurity 
(lower plot) for two different values of ao, as indicated in the figure. The symbols are from MD 
simulations, and the solid (dashed) line the theoretical prediction for ao = 0.95 (0.7) discussed in 
the text. The parameter of the gas are: a = 0.95, N = 359, and W = 50a. The other parameters 
of the impurity are mo = m/2 and ao = a. 

in each of the figures results for two different values of the coefficient of restitution ao are 
reported. 

For the cases reported in Fig. [7| and |Hl the agreement between theory and simulation 
can be considered as quite satisfactory, given that no fitting parameter is being used when 
constructing the density profile for the impurity from the hydrodynamic profiles of the gas. 
On the other hand, for mo = 3m, Fig. |U] shows strong discrepancies between theory and 
simulation results. There are several possible causes for it. One is that the heavy intruder is 
close, on the average, to the vibrating wall and in this region the hydrodynamic description 
for it is not accurate, as indicated in connection with Fig. El In fact, the theoretical prediction 
for Pq{z) is maximum at the wall. Another closely related possible origin of the discrepancy, 
is that we are assuming that the intruder does not affect the hydrodynamic profiles of the 
gas. This may be a bad approximation when the intruder is too massive as compared with 
the gas particles, specially taking into account that the number of the latter is not very 



21 



2 

Nx 



1 



Po 

0.01 

0.005 


100 200 

z/a 

FIG. 8: The same as in Fig. with the only difference that here it is mo = m and the two values 
of ao considered are 0.99 and 0.7, as indicated. 

large. We remind that this number can not be increased too much, in order to keep the 
system fluidized with low density and avoiding the transversal instability mentioned in Sec. 

EH 

We have investigated also the density profile of the intruder for other values of the diam- 
eter ratio o"o/cr, namely 1/2 and 2, and obtained similar results. 

^From Pq(z), the average height of the intruder z can be computed through 

POO 

z = / dzzP (z). (83) 
Jo 

The ratio z /z cm , where z cm is the height of the center of mass of the gas, provides a 
measurement of the relative position of the intruder with respect to the gas. This ratio is 
plotted in Fig. [TU] as a function of for several systems with a = 0.95 and <J /a = 1/2 
(j3 ~ 0.0230). The solid line is the theoretical prediction and the symbols simulation results 
for different values of the mass ratio and «o- Figures ^2 and show the same function but 
for a /a = 1 (0 ~ 0.0172) and a /a = 2 (/3 ~ 0.0115). 

The agreement between theory and simulation is good in all cases, specially in the region 
close to = 1. In fact, for m /m = 3 it is better than it could be expected from Fig. 

This good agreement is not surprising, since we are looking now to the simplest global 
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FIG. 9: The same as in Fig. |HJ with the only difference that here it is ttiq = 3m. 
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FIG. 10: Ratio between the impurity and the gas center of mass positions vs the ratio of the mean 
square velocities </> for systems with gq = cr/2. The solid line is the theoretical prediction discussed 
in the text and the symbols MD simulation results. The same symbol is used for simulations 
differing only on the value of the restitution coefficient ao, which has been varied between 0.3 and 
0.95. The values of the other parameters are given in the text. 
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FIG. 11: The same as in Fig. I1U1 with the only difference that here it is ctq = a. 
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FIG. 12: The same as in Fig. 111! with the only difference that here it is <jq = 2a. 



property of the distribution P${z). On the other hand, the simulation results for cr = 2cr 
are systematically above the theoretical prediction. In general, it seems that the slope of the 
simulation results is smaller than the theoretical prediction, although we have no explanation 
for it. A relevant conclusion to be reached from the above results is the dominant role played 
by the value of in determining the relative position of the intruder with respect to the 
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gas. In an elastic system, is equal to the mass ratio m/mo, but in a granular gas can be 
larger or smaller than the mass ratio due to the lack of energy equipartition. Consequently, 
the relative position of an intruder in a granular gas depends in a rather complicated way 
on the coefficients of restitution of the system, the mass ratio, and the size ratio. 



V. CONCLUSIONS 



In this work, we have studied the hydrodynamic profiles of an impurity immersed in a low 
density granular gas that is in an arbitrary state. Analytical predictions for the temperature 
and density profiles have been derived and compared with molecular dynamics simulations 
for the particular case of a vibrated system. A satisfactory qualitative and quantitative 
agreement has been found, although the discrepancies increase as the difference between the 
masses of the gas particles and the intruder increases. When this is the case, the boundary 
layer next to the vibrating wall becomes rather wide and, moreover, the intruder influences 
the density profile of the gas, contrary to the tracer limit assumed in the theory. 

The theoretical study is based on the validity of the usual Chapman-Enskog procedure to 
generate the hydrodynamic equations of the system when starting from the kinetic (Boltz- 
mann) equations. Then, although the results are quite general and can be applied to a 
variety of situations as the specific one considered here, some caution is needed for other 
particular cases. For instance, in the so-called steady simple shear state, there is an intrinsic 



coupling between inelasticity and gradients 
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271 ] . This implies that the shear rate and 



the coefficient of restitution can not be considered as independent quantities, and the usual 
Chapman-Enskog procedure becomes inaccurate. This is in contrast with the situation con- 
sidered here, in which the gradients are controlled by the vibrating wall and the external 
gravitational field. In general, each physical state should be checked to verify whether a 
direct expansion in the gradients is legitimate. 

For the ratio between the local temperatures of the intruder and the gas, the theory 
predicts that it is constant in the bulk of the system. This has been confirmed by the MD 
simulations. The same behavior has been found in experiments j?| Q and MD simulations 
of vibrated systems, outside the dilute and tracer limits. Nevertheless, an extension of 

28| . indicates that the ratio depends also on the local 



the analysis presented in this paper 
densities ratio. What happens is that this dependence is rather weak and gives corrections 
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inside the statistical uncertainties for the values of the parameters considered in the above 
mentioned works. 

A qualitative discussion of the relative position of the intruder with respect to the gas was 
given in [12j , where the crucial role played by the non-equipartition of energy was discussed. 
More precisely, it was shown that the intruder was above or below the center of mass of 
the gas depending on the sign of the pressure diffusion coefficient D p , that in turn is mostly 
determined by the sign of <p — 1, as seen in Eq. (|75p. This is in agreement with the simulation 
results reported in Figs. ITUHT21 It can be concluded that the difference in temperature of the 
components of a mixture must be taken into account when studying segregation phenomena, 
i.e. the demix of the components when shaken. Of course, at higher densities and finite 
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301 ] are also important and 



volume fraction of the components, other mechanisms 
even dominant. 

It is worth to mention that the work reported here provides another example of the ac- 
curacy of hydrodynamics to describe granular systems under rather extreme conditions of 
inelasticity. Although the inelasticity of the gas has been kept small in order to avoid the 
transversal instability, predicted by the hydrodynamic theory itself, the coefficient of inelas- 
ticity for the gas-impurity collisions has been varied over a wide range (typically between 0.3 
and 0.99) and a reasonable good agreement between theory and simulation has been found. 
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APPENDIX: HYDRODYNAMIC PROFILES OF THE GAS 

In ref. 



2l| . analytical expressions for the density and temperature profiles of an open 
vibrated granular gas in the steady state were derived. They are given by 

C d a d ^Ja(a) [AI,(0 + BK^)Y 



n \0 = - - ; r— — ~ — ~ l2 ' ( A - x ) 



T i / 2 (o = [A/ 7 (o + ^ 7 (0]r 7 , (a.2) 
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where = 2 a/2 for d = 2 and Cd = tt\^2 for d = 3, A and B are constants to be determined 
from some boundary conditions, and I 7 and denote the modified Bessel functions of first 
and second kind, respectively. The dimensionless length scale £ is defined by 

f =^r^- <A - 3) 

with X(z) = (Cdna d ~ 1 )~ 1 being the local mean free path. Finally, a(a) and 7(a) are functions 
of the coefficient of normal restitution, both vanishing in the elastic limit a — > 1. Their 



explicit expressions are given in 



2l\ . Note that the maximum value of £ is 



6 = ^pv, (A . 4) 



N being the total number of particles and W the length or area of the vibrating wall. 
Using the hydrodynamic relation 

dp 



-nmg , (A.5) 



we can write g(z), defined in Eq. (jSTjl . as 

^) = -r("^ + W&j- (A - 6) 

The expression for dT/dz is easily obtained by using Eqs. (|A.l|) - (jA.3J) and the properties of 



the modified Bessel functions 



3l|, 



dT = 2mg {j [A/ 7+1 (0 - BK 1+1 {j)\ 

dz AI^+BK^i) • 1 • ' 

In this way, we have an explicit expression for the integral appearing in Eq. ()82|) that, 
therefore, can be evaluated numerically for each set of values of the parameters. Of course, 
this requires to determine first the constants A and B. In the results to be reported here, 
it has been done from the simulation data themselves, namely by fitting the position and 
value of the temperature minimum 2l|. Another possibility that has been proven to give 
accurate results would be to require the hydrodynamic heat flux to vanish at infinite height 
j^ . Then for a granular gas with a = 0.95, = 359, Vw = 5, and W = 50, leading to 
7 ~ 0.021, a ~ 0.010, and £ - 2.065, the values A ~ 5.41 and B ~ 0.14 are found. The 
solid lines in Fig. ^ are the fit to the hydrodynamic profiles of the gas obtained in this way. 
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Once A and B are determined, it is possible to evaluate numerically the right hand side 
of Eq. ()82j) for each set of values of the parameters defining the mechanical properties of the 
intruder. 
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